Projektuppgift

Linjär Algebra

Högskolan i Halmstad

Dawid Ejdeholm, Jacob Harsten, Harald Lilja & Jonatan Jönsson

Uppgift 1

Clear["Global`*"]
(OverVector[(Subscript[e, x])])=UnitVector[3,1];
(OverVector[(Subscript[e, y])]) = UnitVector[3,2];
(OverVector[(Subscript[e, z])])=UnitVector[3,3];
m = 150;
k = 2000;
OG = 1.4;
OB = 2.1;
OC = 4.2;
d = (({{6, 2}}));
g = 9.8;
"Datan laddad"

a)  Förklara varför det här är lämpligt att beräkna momenten med avseende på punkten O.

Eftersom hävarmen är fäst i punkten O.

b)  Identifiera alla krafter som ger ett moment (m.a.p. O) och uttryck både krafterna och ortsvektorer för
krafternas angrepspunkter på komponentform.

Tyngdkraft - mg, punkt G
Fjäderkraft - ks, punkt B
Spännkraft - T, punkt C
Ortsvektorerna:

Clear[θ];
θ = 90;
(OverVector[OG])=1.4*Sin[θ °]*(OverVector[(Subscript[e, x])])-1.4*Cos[θ °]*(OverVector[(Subscript[e, y])])+0*(OverVector[(Subscript[e, z])]);
(OverVector[OB])=2.1*Sin[θ °]*(OverVector[(Subscript[e, x])])-2.1*Cos[θ °]*(OverVector[(Subscript[e, y])])+0*(OverVector[(Subscript[e, z])]);
(OverVector[OC])=4.2*Sin[θ °]*(OverVector[(Subscript[e, x])])-4.2*Cos[θ °]*(OverVector[(Subscript[e, y])])+0*(OverVector[(Subscript[e, z])]);
"Ortsvektorerna laddad"

Tyngdkraft:

(OverVector[(Subscript[F, mg])])=0*(OverVector[(Subscript[e, x])])-m*g*(OverVector[(Subscript[e, y])])+0*(OverVector[(Subscript[e, z])])

Kryssprodukten ger momentet med rätt riktning, i detta fall negativ z-axel enligt högerhandsregeln.

(OverVector[(Subscript[M, mg])])=Cross[(OverVector[OG]),(OverVector[(Subscript[F, mg])])]

Fjäderkraft:
	Kraftvektor = Belopp av kraft × normerad riktningsvektor
	(OverVector[(Subscript[F, fj])])=k*s*(OverVector[(Subscript[n, fj])])
		(OverVector[(Subscript[n, fj])])=((OverVector[BA])/(Norm[(OverVector[BA])]))
			(OverVector[BA])=(OverVector[OA])-(OverVector[OB])

(OverVector[OA])=0*(OverVector[(Subscript[e, x])])-4.2*(OverVector[(Subscript[e, y])])+0*(OverVector[(Subscript[e, z])]);
(OverVector[BA])=(OverVector[OA])-(OverVector[OB]);
(OverVector[(Subscript[n, fj])])= ((OverVector[BA])/(Norm[(OverVector[BA])]));
s=Norm[(OverVector[BA])]-(OC-OB);
(OverVector[(Subscript[F, fj])])= k*s*(OverVector[(Subscript[n, fj])])

Momentet:

(OverVector[(Subscript[M, fj])])=Cross[(OverVector[OB]),(OverVector[(Subscript[F, fj])])]

Spännkraft:
	(OverVector[(Subscript[F, sp])])=T*(OverVector[(Subscript[n, sp])])
		(OverVector[(Subscript[n, sp])])=((OverVector[CD])/(Norm[(OverVector[CD])]))
			(OverVector[CD])=(OverVector[OD])-(OverVector[OC])

(OverVector[OD])=6*(OverVector[(Subscript[e, x])])+2*(OverVector[(Subscript[e, y])])+0*(OverVector[(Subscript[e, z])]);
(OverVector[CD])=(OverVector[OD])-(OverVector[OC]);
(OverVector[(Subscript[n, sp])])=((OverVector[CD])/(Norm[(OverVector[CD])]));

Vi vet att momentet som spännkraften ger är lika stort som (OverVector[(Subscript[M, fj])])+(OverVector[(Subscript[M, mg])])
Vi kan då lösa ut T

(OverVector[(Subscript[M, sp])])=Norm[(OverVector[(Subscript[M, fj])])+(OverVector[(Subscript[M, mg])])];
T = ((OverVector[(Subscript[M, sp])])/(Norm[Cross[(OverVector[(Subscript[n, sp])]),(OverVector[OC])]]))

Kontroll:

(OverVector[(Subscript[M, sp])]);
Simplify[Cross[T*(OverVector[(Subscript[n, sp])]),(OverVector[OC])]]

c) Ställ upp och lös jämviktsekvationen med avseende på spännkraften i linan, T.
Uttryck denna som en funktion av vinkeln T(θ).

Ifrån b) har vi redan:

Clear[T];
T[θ_]:=((Sqrt[0.` +((Abs[0.` -19698.` Sin[° θ]+((37044.` Sin[° θ])/(Sqrt[((Abs[-4.2`+2.1` Cos[° θ]])^2)+((Abs[0.` -2.1` Sin[° θ]])^2)]))])^2)])/(Sqrt[0.` +((Abs[0.` -((25.200000000000003` Cos[° θ])/(Sqrt[((Abs[2.` +4.2` Cos[° θ]])^2)+((Abs[6.` -4.2` Sin[° θ]])^2)]))-((8.4` Sin[° θ])/(Sqrt[((Abs[2.` +4.2` Cos[° θ]])^2)+((Abs[6.` -4.2` Sin[° θ]])^2)]))])^2)]))
T[45]

d) Visa T(θ) grafiskt.

Vi plottar T[θ] med (Subscript[θ, min])=0 och (Subscript[θ, max])=90°+ArcTan[2/6]
Vi använder detta (Subscript[θ, max])eftersom det är så långt hävarmen kan gå genom att dra i linan.

Plot[((Sqrt[0.` +((Abs[0.` -19698.` Sin[θ °]+((37044.` Sin[θ °])/(Sqrt[((Abs[-4.2`+2.1` Cos[θ °]])^2)+((Abs[0.` -2.1` Sin[θ °]])^2)]))])^2)])/(Sqrt[0.` +((Abs[0.` -((25.200000000000003` Cos[θ °])/(Sqrt[((Abs[2.` +4.2` Cos[θ °]])^2)+((Abs[6.` -4.2` Sin[θ °]])^2)]))-((8.4` Sin[θ °])/(Sqrt[((Abs[2.` +4.2` Cos[θ °]])^2)+((Abs[6.` -4.2` Sin[θ °]])^2)]))])^2)])),{θ,0,N[90+((ArcTan[(2/6)])/°)]}]

e) Antag att linan maximalt tål en belastning T = 5.0 kN.
Vid vilken vinkel når vi detta värde för spännkraften ?

Vi använder våran funktion T[θ]

Clear[θ];
Solve[T[θ]==5000&& θ>0&&θ<=N[90+((ArcTan[(2/6)])/°)],θ]

Vid θ = 95.8198°

f) Vad blir det maximala värdet för θ om vi inte behöver beakta linans begränsningar ?

Detta är samma vinkel som vi tidigare använda som (Subscript[θ, max ]), alltså:

(Subscript[θ, max])=N[90+((ArcTan[(2/6)])/°)]

Uppgift 2

Ett experiment i en accelerator producerar ett grund¨amne som enbart finns i form av tre kortlivade isotoper.
Dessa har halveringstider 2, 7 respektive 24 minuter. Den totala mängden/massan av ämnet, M, registereras
sedan under någon tid efter experimentet. Se tabellen nedan.

Remove["Global`*"]

a)  Använd minsta kvadratmetoden och bestäm approximativa värden för ursprunglig mängd för var och en
av isotoperna.

M(t) = A(e^(((-ln 2)/2)t))+B(e^(((-ln 2)/7)t))+C(e^(((-ln 2)/24)t))

𝔸 = (({{(2^((-1)/2)), (2^((-1)/7)), (2^((-1)/24))}, {(2^((-2)/2)), (2^((-2)/7)), (2^((-2)/24))}, {(2^((-3)/2)), (2^((-3)/7)), (2^((-3)/24))}, {(2^((-4)/2)), (2^((-4)/7)), (2^((-4)/24))}, {(2^((-5)/2)), (2^((-5)/7)), (2^((-5)/24))}, {(2^((-6)/2)), (2^((-6)/7)), (2^((-6)/24))}, {(2^((-7)/2)), (2^((-7)/7)), (2^((-7)/24))}, {(2^((-8)/2)), (2^((-8)/7)), (2^((-8)/24))}, {(2^((-9)/2)), (2^((-9)/7)), (2^((-9)/24))}, {(2^((-10)/2)), (2^((-10)/7)), (2^((-10)/24))}}));
𝕐 = (({{12.57}, {10.22}, {8.42}, {6.92}, {5.96}, {5.09}, {5.43}, {4.01}, {3.62}, {3.28}}));
𝕏 = (({{a}, {b}, {c}}));
sol = Solve[𝔸\[Transpose].𝔸.𝕏== 𝔸\[Transpose].𝕐]

data = {{1,12.57}, {2,10.22},{3, 8.42},{4,6.92},{5,5.96},{6,5.09},{7,5.43},{8,4.01},{9,3.62},{10,3.28}};
x = (2^((-t)/2));
y = (2^((-t)/7));
z= (2^((-t)/24));
FindFit[data,a*x+b*y+c*z,{a,b,c},t]

Svar: 
a-> 9.2103
b-> 4.7243
c-> 1.80546

b)  De framräknade värdena ovan ger en modell, Ma(t), för totala mängden som funktion av tiden.
Gör en grafisk jämförelse av Ma(t) och M i intervallet 0 ≤ t ≤ 10.

Remove["Global`*"]

Plot[9.2103*(2^((-t)/2))+4.7243*(2^((-t)/7))+1.80546*(2^((-t)/24)),{t,0,10},Epilog->{PointSize[.012],
Red,Point[data]},AxesLabel->{"t","M"},
LabelStyle->Directive[Italic,13]]

c) Antag i stället att två av isotoperna har samma halveringstid. Förklara varför det nu blir omöjligt att
använda minsta kvadratmetoden (om vi inte har ytterligare information om initiala förhållandet mellan
mängderna av de två isotoper som har samma halveringstid).

𝔸 = (({{(2^((-1)/2)), (2^((-1)/2)), (2^((-1)/24))}, {(2^((-2)/2)), (2^((-2)/2)), (2^((-2)/24))}, {(2^((-3)/2)), (2^((-3)/2)), (2^((-3)/24))}, {(2^((-4)/2)), (2^((-4)/2)), (2^((-4)/24))}, {(2^((-5)/2)), (2^((-5)/2)), (2^((-5)/24))}, {(2^((-6)/2)), (2^((-6)/2)), (2^((-6)/24))}, {(2^((-7)/2)), (2^((-7)/2)), (2^((-7)/24))}, {(2^((-8)/2)), (2^((-8)/2)), (2^((-8)/24))}, {(2^((-9)/2)), (2^((-9)/2)), (2^((-9)/24))}, {(2^((-10)/2)), (2^((-10)/2)), (2^((-10)/24))}}));
𝕐 = (({{12.57}, {10.22}, {8.42}, {6.92}, {5.96}, {5.09}, {5.43}, {4.01}, {3.62}, {3.28}}));
𝕏 = (({{a}, {b}, {c}}));
sol = Solve[𝔸\[Transpose].𝔸.𝕏== 𝔸\[Transpose].𝕐]

Uppgift 3

Remove["Global`*"];

a) Förklara varför (Subscript[n, 0])(t) = (Subscript[f, 3])(Subscript[n, 3])(t) + (Subscript[f, 4])(Subscript[n, 4])(t) + (Subscript[f, 5])(Subscript[n, 5])(t) och (Subscript[n, k+1])(t + 1) = (Subscript[s, k])(Subscript[n, k])(t), k = 0, 1, 2, 3, 4, där t ≥ 0.

(Subscript[n, 0])(t) = (Subscript[f, 3])(Subscript[n, 3])(t) + (Subscript[f, 4])(Subscript[n, 4])(t) + (Subscript[f, 5])(Subscript[n, 5])(t)
Antalet nyfödda = Antalet honor som föds per individ * Antalet honor 

(Subscript[n, k+1])(t+1)=(Subscript[s, k])(Subscript[n, k])(t),k=0,1,2,3,4,där t≥0.
Nästkommande års population = Överlevnadsandel * Antalet individer

b)  Formulera resultaten ovan som en 5 × 5 matris P verkande på en kolonnmatris

ℕ=(({{(Subscript[n, 1])}, {(Subscript[n, 2])}, {(Subscript[n, 3])}, {(Subscript[n, 4])}, {(Subscript[n, 5])}}));
matrisen=(({{0, 0, (Subscript[f, 3]), (Subscript[f, 4]), (Subscript[f, 5])}, {(Subscript[s, 1]), 0, 0, 0, 0}, {0, (Subscript[s, 2]), 0, 0, 0}, {0, 0, (Subscript[s, 3]), 0, 0}, {0, 0, 0, (Subscript[s, 4]), 0}}));
N(t+1)MatrixForm[matrisen.ℕ]

c) Bestäm populationen vid tiden t = 30 (år). Går det att se en tendens ?

ℕ=(({{5}, {2}, {2}, {1}, {1}}))(10^3);
P=(({{0, 0, 3/2, 4/2, 3/2}, {2/5, 0, 0, 0, 0}, {0, 3/5, 0, 0, 0}, {0, 0, 7/10, 0, 0}, {0, 0, 0, 1/2, 0}}));
antal =MatrixForm[P.ℕ]

λ=0.94879;
egenvektor = {{9.6472},{4.06716},{2.57201},{1.89758},{1}};
P.egenvektor
λ*egenvektor
"P * egenvektorn = lambda * egenvektor";

MatrixForm[N[MatrixPower[P,30].ℕ]]

d) Visa att populationen kommer dö ut.

ℕ=(({{5}, {2}, {2}, {1}, {1}}))(10^3);
P=(({{0, 0, 3/2, 4/2, 3/2}, {2/5, 0, 0, 0, 0}, {0, 3/5, 0, 0, 0}, {0, 0, 7/10, 0, 0}, {0, 0, 0, 1/2, 0}}));
Limit[MatrixPower[P,x].ℕ, x->∞]

e) Antag att det är praktiskt genomförbart att lyfta överlevnadsandelen för årsklass 0 med t.ex. riktad
foderutsättning. Vad måste (Subscript[s, 0]) ändras till för att få en stabil population ?

P=(({{0, 0, S*3, S*4, S*3}, {2/5, 0, 0, 0, 0}, {0, 3/5, 0, 0, 0}, {0, 0, 7/10, 0, 0}, {0, 0, 0, 1/2, 0}}));
Solve[Det[IdentityMatrix[5]-P]==0,S]

f) Välj det nya värdet för s0 från deluppgiften ovan och antag att vi startar med samma ursprungliga population
som ovan. Vad blir totala antalet honor efter lång tid (inkluderad åldersklass 0) ?

ℕ=(({{5}, {2}, {2}, {1}, {1}}))(10^3);
S=(250/411);
P=(({{0, 0, S*3, S*4, S*3}, {2/5, 0, 0, 0, 0}, {0, 3/5, 0, 0, 0}, {0, 0, 7/10, 0, 0}, {0, 0, 0, 1/2, 0}}));

Total[N[MatrixPower[P,43].ℕ]]

Svar: Efter 43 år finns det totalt 11969 st honor.